R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

library(raster)
size_dat=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/fire_growth_5days_v4.txt", header=T, row.names=NULL)

size_dat=as.data.frame(size_dat)
size_dat=size_dat[-1,]
size_dat=size_dat[,-1]
colnames(size_dat)=c("firename","year","cause","size1","size2","size3","size4","size5","final_firesize","mean_precip1","mean_precip2","mean_precip3","mean_precip4","mean_precip5","mean_tmax1","mean_tmax2","mean_tmax3","mean_tmax4","mean_tmax5","mean_tmean1","mean_tmean2","mean_tmean3","mean_tmean4","mean_tmean5","mean_vpdmax1","mean_vpdmax2","mean_vpdmax3","mean_vpdmax4","mean_vpdmax5","mean_windspeed1","mean_windspeed2","mean_windspeed3","mean_windspeed4","mean_windspeed5","landcover","ecosystem","biomass","elevation")

size_dat2 <- data.frame(lapply(size_dat, function(x) as.numeric(as.character(x))))

NAs introduced by coercion

size_dat2$human = 0
size_dat2$human[size_dat2$cause !=1 & size_dat2$cause !=14 & size_dat2$cause !=17]=1
size_dat2$human[size_dat2$cause ==1 ]=2

pro1 =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1 ), ]
pro2 =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1 ), ]

t.test(pro1$size1,pro2$size1)

    Welch Two Sample t-test

data:  pro1$size1 and pro2$size1
t = -2.0397, df = 74.381, p-value = 0.04493
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.0986230  -0.2593373
sample estimates:
mean of x mean of y 
 3.056385 14.235365 

t.test(pro1$size2,pro2$size2)

    Welch Two Sample t-test

data:  pro1$size2 and pro2$size2
t = -2.758, df = 73.138, p-value = 0.007341
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -32.702045  -5.266208
sample estimates:
mean of x mean of y 
 8.254394 27.238521 

t.test(pro1$size3,pro2$size3)

    Welch Two Sample t-test

data:  pro1$size3 and pro2$size3
t = -2.9002, df = 58.203, p-value = 0.005254
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -51.965993  -9.526741
sample estimates:
mean of x mean of y 
 14.10484  44.85121 

t.test(pro1$size4,pro2$size4)

    Welch Two Sample t-test

data:  pro1$size4 and pro2$size4
t = -3.1797, df = 53.078, p-value = 0.002461
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -74.99798 -16.98040
sample estimates:
mean of x mean of y 
 17.68430  63.67349 

t.test(pro1$size5,pro2$size5)

    Welch Two Sample t-test

data:  pro1$size5 and pro2$size5
t = -3.4218, df = 39.443, p-value = 0.001462
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -117.71273  -30.26844
sample estimates:
mean of x mean of y 
 19.71795  93.70853 

pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1), ]
length(pro$year)

[1] 68

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1), ]
length(pro$year)

[1] 77

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 2), ]
length(pro$year)

[1] 18

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 2), ]
length(pro$year)

[1] 9

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

length(!is.na(pro$size5))

[1] 9

pro =size_dat2[which(size_dat2$human == 1 & size_dat2$ecosystem == 6), ]
length(pro$year)

[1] 41

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 2& size_dat2$ecosystem == 6), ]
length(pro$year)

[1] 79

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

NA
NA

plot fire size map QGIS

library(data.table)
DT= data.table(res)
fire_size = DT[ , .SD[which.min(growth_km)], by = firename]

fire_size = fire_size[,c()]
spread_list =list.files(viirs_dir, pattern = "_daily.shp$", recursive = TRUE, full.names=T)
i=1
 l=shapefile(spread_list[i])
for (i in 2:length(spread_list)){
  p=shapefile(spread_list[i])
  l <- rbind(l, p, makeUniqueIDs = TRUE) 
}
 l$human = 0
l$human[l$cause !=1 & l$cause !=14 & l$cause !=17]=1
l$human[l$cause ==1 ]=2

     writeOGR(l, "/Users/stijnhantson/Documents/projects/VIIRS_ros/", layer= "all_fires", driver="ESRI Shapefile", overwrite_layer = T)

extract number and size statistics from frap

library(raster)

Loading required package: sp

dr =shapefile("/Users/stijnhantson/Documents/projects/VIIRS_ros/frap_subset.shp")
dr1 =shapefile("/Users/stijnhantson/Documents/data/FRAP/fire18_1.shp")
dr1$YEAR_=as.numeric(as.character(dr1$YEAR_))
dr1$Shape_Area=as.numeric(as.character(dr1$Shape_Area))
dr1=dr1[!is.na(dr1$YEAR_), ]
dr1=dr1[dr1$YEAR_>2011,]
sum(dr1$Shape_Area)

[1] 25651974971

sum(dr$Shape_Area)

[1] 22801510990

1. frap
2. only fire growth datasets

• prepare final dataset to open

daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/final_dataset_V3.txt",header=T)

res=as.data.frame(daily_res)

res$mean_ros =as.numeric(as.character(res$mean_ros))
res$max_ros =as.numeric(as.character(res$max_ros))
res$median95_ros =as.numeric(as.character(res$median95_ros))
res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
#res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$growth =as.numeric(as.character(res$growth))
res$total_area =as.numeric(as.character(res$total_area))
res$mean_frp =as.numeric(as.character(res$mean_frp))
res$frp_95 =as.numeric(as.character(res$frp_95))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))
res$biomass =as.numeric(as.character(res$biomass))
res$year =as.numeric(as.character(res$year))
res$month =as.numeric(as.character(res$month))
res$doy_out =as.numeric(as.character(res$doy_out))

res = res[-1,]
res$per_ba = res$growth/res$total_area
res$growth_km =res$growth/1000000

res$human = 0
res$human[res$cause !=1 & res$cause !=14 & res$cause !=17]=1
res$human[res$cause ==1 ]=2

res$ros_km = (res$median95_ros*24)/1000
res$ros_mean_km = (res$mean_ros*24)/1000

plot mean vs max fire rate-of-spread

#summary(res)
plot(res$ros_km,res$ros_mean_km, xlab="maximum fire-spread-rate (km/day",ylab="mean fire-spread-rate (km/day)")


tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/mean_vs_max_ros_v1.tif", width = 5, height = 5, units = 'in', res = 300)

plot(res$ros_km,res$ros_mean_km, xlim=c(0,25),ylim=c(0,10), xlab="maximum fire rate-of-spread (km/day)",ylab="mean fire rate-of-spread (km/day)", cex.lab=1.3,cex.axis = 1.25)
dev.off()

quartz_off_screen 
                2 

difference between human and lightnign fires

me=0
me1=0
days = c("day1","day2","day3","day4","day5")
days1 = c("1","2","3","4","5")
pro1 = res[res$fire_day == 1 & res$human == 1,]
pro2 = res[res$fire_day == 2 & res$human == 1,]
pro3 = res[res$fire_day == 3 & res$human == 1,]
pro4 = res[res$fire_day == 4 & res$human == 1,]
pro5 = res[res$fire_day == 5 & res$human == 1,]
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

pro1h = res[res$fire_day == 1 & res$human == 2,]
pro2h = res[res$fire_day == 2 & res$human == 2,]
pro3h = res[res$fire_day == 3 & res$human == 2,]
pro4h = res[res$fire_day == 4 & res$human == 2,]
pro5h = res[res$fire_day == 5 & res$human == 2,]
me1[1] =mean(pro1h$growth_km,na.omit=T)
me1[2] =mean(pro2h$growth_km,na.omit=T)
me1[3] =mean(pro3h$growth_km,na.omit=T)
me1[4] =mean(pro4h$growth_km,na.omit=T)
me1[5] =mean(pro5h$growth_km,na.omit=T)

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure1_v1.tif", width = 10, height = 5, units = 'in', res = 300)

par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^-1*')'),ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)

dev.off()

quartz_off_screen 
                2 

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_figure1_v1.tif", width = 10, height = 5, units = 'in', res = 300)

par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^-1*')'),ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)

dev.off() 

quartz_off_screen 
                2 

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/all.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

t.test(log(pro1$growth_km),log(pro1h$growth_km))

    Welch Two Sample t-test

data:  log(pro1$growth_km) and log(pro1h$growth_km)
t = 6.5083, df = 169.36, p-value = 8.265e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.492386 2.791866
sample estimates:
 mean of x  mean of y 
 1.3226886 -0.8194377 

t.test(log(pro2$growth_km),log(pro2h$growth_km))

    Welch Two Sample t-test

data:  log(pro2$growth_km) and log(pro2h$growth_km)
t = 6.9843, df = 134.4, p-value = 1.199e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.369607 2.451706
sample estimates:
mean of x mean of y 
2.7559564 0.8453003 

t.test(log(pro3$growth_km),log(pro3h$growth_km))

    Welch Two Sample t-test

data:  log(pro3$growth_km) and log(pro3h$growth_km)
t = 5.3694, df = 137.17, p-value = 3.286e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.9293526 2.0129234
sample estimates:
mean of x mean of y 
 3.186372  1.715234 

t.test(log(pro4$growth_km),log(pro4h$growth_km))

    Welch Two Sample t-test

data:  log(pro4$growth_km) and log(pro4h$growth_km)
t = 4.6217, df = 122.92, p-value = 9.467e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.6909157 1.7261008
sample estimates:
mean of x mean of y 
 3.513204  2.304696 

t.test(log(pro5$growth_km),log(pro5h$growth_km))

    Welch Two Sample t-test

data:  log(pro5$growth_km) and log(pro5h$growth_km)
t = 4.9328, df = 103.98, p-value = 3.089e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.9170094 2.1499736
sample estimates:
mean of x mean of y 
 3.895266  2.361774 

##################3 for western cordillera ecoregion ##################

me=0
me1=0
pro1 = res[res$fire_day == 1 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

t.test(log(pro1$growth_km),log(pro1h$growth_km))

    Welch Two Sample t-test

data:  log(pro1$growth_km) and log(pro1h$growth_km)
t = 3.2086, df = 69.426, p-value = 0.002019
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.5363985 2.2992588
sample estimates:
 mean of x  mean of y 
 0.4219047 -0.9959240 

t.test(log(pro2$growth_km),log(pro2h$growth_km))

    Welch Two Sample t-test

data:  log(pro2$growth_km) and log(pro2h$growth_km)
t = 4.9267, df = 97.316, p-value = 3.427e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.9938201 2.3346153
sample estimates:
mean of x mean of y 
 2.303104  0.638886 

t.test(log(pro3$growth_km),log(pro3h$growth_km))

    Welch Two Sample t-test

data:  log(pro3$growth_km) and log(pro3h$growth_km)
t = 3.4468, df = 82.843, p-value = 0.0008936
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4842092 1.8055445
sample estimates:
mean of x mean of y 
 2.695279  1.550402 

t.test(log(pro4$growth_km),log(pro4h$growth_km))

    Welch Two Sample t-test

data:  log(pro4$growth_km) and log(pro4h$growth_km)
t = 3.019, df = 73.938, p-value = 0.003479
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.3104149 1.5156656
sample estimates:
mean of x mean of y 
 3.092997  2.179957 

t.test(log(pro5$growth_km),log(pro5h$growth_km))

    Welch Two Sample t-test

data:  log(pro5$growth_km) and log(pro5h$growth_km)
t = 3.1504, df = 55.57, p-value = 0.002625
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4138127 1.8597722
sample estimates:
mean of x mean of y 
 3.402659  2.265867 

mean(pro1$growth_km,na.omit=T)

[1] 16.60532

mean(pro1h$growth_km,na.rm=T)

[1] 2.714784

for mediteranean california

pro1 = res[res$fire_day == 1 & res$human == 1 & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & (res$eco1 == 11),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & (res$eco1 == 11),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

for difference in autumn

pro1 = res[res$fire_day == 1 & res$human == 1 & res$month >8,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month >8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month >8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month >8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month >8,]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & res$month >8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month >8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month >8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month >8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month >8,]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/autumn.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,250),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

for difference in summer

pro1 = res[res$fire_day == 1 & res$human == 1 & res$month <=8 ,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month <=8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month <=8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month <=8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month <=8,]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & res$month <=8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month <=8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month <=8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month <=8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month <=8,]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/summer.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

for difference in summer in western cordillera

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

t.test(log(pro1$growth_km),log(pro1h$growth_km))

    Welch Two Sample t-test

data:  log(pro1$growth_km) and log(pro1h$growth_km)
t = 1.3062, df = 51.551, p-value = 0.1973
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.3236595  1.5299715
sample estimates:
 mean of x  mean of y 
-0.1773847 -0.7805406 

t.test(log(pro2$growth_km),log(pro2h$growth_km))

    Welch Two Sample t-test

data:  log(pro2$growth_km) and log(pro2h$growth_km)
t = 3.6535, df = 78.926, p-value = 0.0004639
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.5639969 1.9140943
sample estimates:
mean of x mean of y 
2.0787175 0.8396718 

t.test(log(pro3$growth_km),log(pro3h$growth_km))

    Welch Two Sample t-test

data:  log(pro3$growth_km) and log(pro3h$growth_km)
t = 2.1785, df = 57.055, p-value = 0.03352
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.05727942 1.36024900
sample estimates:
mean of x mean of y 
 2.551723  1.842958 

t.test(log(pro4$growth_km),log(pro4h$growth_km))

    Welch Two Sample t-test

data:  log(pro4$growth_km) and log(pro4h$growth_km)
t = 2.353, df = 57.21, p-value = 0.02208
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1074894 1.3347908
sample estimates:
mean of x mean of y 
 2.992332  2.271192 

t.test(log(pro5$growth_km),log(pro5h$growth_km))

    Welch Two Sample t-test

data:  log(pro5$growth_km) and log(pro5h$growth_km)
t = 2.3446, df = 41.508, p-value = 0.02391
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.105315 1.410281
sample estimates:
mean of x mean of y 
 3.313679  2.555881 

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_summer.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

for difference in autumn in western cordillera

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_autumn.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,250),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

for difference in summer in meditereanean

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

#t.test(log(pro1$growth_km),log(pro1h$growth_km))
#t.test(log(pro2$growth_km),log(pro2h$growth_km))
#t.test(log(pro3$growth_km),log(pro3h$growth_km))
#t.test(log(pro4$growth_km),log(pro4h$growth_km))
#t.test(log(pro5$growth_km),log(pro5h$growth_km))

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_summer.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

` ################## for difference in autumn in mediteranean ##################

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_autumn.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,400),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                2 

how many days does it take to reach 75% burnt area

res75 = res[res$per_ba > 0.75,]
#peak_day = as.data.frame(aggregate(res75$fire_day, by = list(res75$firename,res75$cause), min))
#peak_day=subset(peak_day,x < 55)
#hi = hist(peak_day$x,prob =F, breaks= c(0:54), xlim=c(0,55), ylab="number of fires", xlab="days", cex.lab=1.4,cex.axis=1.3)

out1 = subset(res75,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 )
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak=as.data.frame(aggregate(res75$fire_day, by = list(res75$firename), min))
quantile(peak_day1$x,0.50,type=3) 

50% 
 10 

quantile(peak_day2$x,0.50,type=3) 

50% 
  3 

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/time_to_reach75.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  mediteranean  ###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11)
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,15))

axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  north cal  ###########
out1 = subset(res75,cause == 1  & (res75$eco1 == 6 | res75$eco1 == 7))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14  & (res75$eco1 == 6 | res75$eco1 == 7))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,15))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  mediteranean  SUMMER ###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11 & ( res75$month <=8 & res75$month >5 ))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11 & ( res75$month <=8 & res75$month >5 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med_summer.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,10))

axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  north cal  summer ###########
out1 = subset(res75,cause == 1  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month <=8 & res75$month >5 ))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month <=8 & res75$month >5 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north_summer.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,12))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  mediteranean  autumn###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11 & ( res75$month >8))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11 & ( res75$month >8 ))
#peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

#peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
#hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(0,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med_autumn.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,5))

axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  north cal  autumn ###########
out1 = subset(res75,cause == 1  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month >8))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month >8 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north_autumn.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,5))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

res=res[res$ros_km>0,]
res_f = res[res$max_land == 1,]
res_p = res[res$max_land != 1,]

#just show the plot here
plot(res_f$mean_frp~res_f$ros_km,log="xy",xlim=c(0.005,30),ylim=c(0.1,180),xaxt="n",ylab="mean FRP (MW)",xlab="Rate-of-Spread (km/day)", cex.lab=1.4,cex.axis = 1.3,col="darkgreen")

marks=c(0.01,0.1,1,10)
marks1=c(0.1,0.5,5,50)
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/fig_FRP_ros_v3.tif", width = 5, height = 5, units = 'in', res = 300)

plot(res_f$mean_frp~res_f$ros_km,log="xy",xlim=c(0.005,30),ylim=c(0.1,180),xaxt="n",yaxt="n",ylab="mean FRP (MW)",xlab="Rate-of-Spread (km/day)", cex.lab=1.4,cex.axis = 1.3,col="darkgreen")
points(res_p$mean_frp~res_p$ros_km,col="orange")
axis(1,at=marks,labels=marks,cex.axis=1.4 )
axis(2,at=marks1,labels=marks1,cex.axis=1.4 )
legend( x="topleft",legend=c("Forest","Grass & shrub"),col=c("darkgreen","orange"),cex=1.2,pch=1,bty = "n")
dev.off()

quartz_off_screen 
                2 

difference in fire size for first 5 days across california and both ecosystems

 
res$ros1 = res$max_ros+1


out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

hum1 = out2[out2$fire_day ==1,]
hum2 = out2[out2$fire_day ==2,]
hum3 = out2[out2$fire_day ==3,]
hum4 = out2[out2$fire_day ==4,]
hum5 = out2[out2$fire_day ==5,]
lig1 = out1[out1$fire_day ==1,]
lig2 = out1[out1$fire_day ==2,]
lig3 = out1[out1$fire_day ==3,]
lig4 = out1[out1$fire_day ==4,]
lig5 = out1[out1$fire_day ==5,]
mean(hum1$growth)

[1] 18753898

mean(hum2$growth)

[1] 36707455

mean(hum3$growth)

[1] 57176147

mean(hum4$growth)

[1] 80694648

mean(hum5$growth)

[1] 115079142

mean(lig1$growth)

[1] 2875395

mean(lig2$growth)

[1] 10065897

mean(lig3$growth)

[1] 18876394

mean(lig4$growth)

[1] 28689523

mean(lig5$growth)

[1] 37533565

t.test(log10(hum1$growth),log10(lig1$growth))

    Welch Two Sample t-test

data:  log10(hum1$growth) and log10(lig1$growth)
t = 6.5083, df = 169.36, p-value = 8.265e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.6481351 1.2124922
sample estimates:
mean of x mean of y 
 6.574436  5.644123 

t.test(log10(hum2$growth),log10(lig2$growth))

    Welch Two Sample t-test

data:  log10(hum2$growth) and log10(lig2$growth)
t = 6.9843, df = 134.4, p-value = 1.199e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.5948126 1.0647622
sample estimates:
mean of x mean of y 
 7.196897  6.367109 

t.test(log10(hum3$growth),log10(lig3$growth))

    Welch Two Sample t-test

data:  log10(hum3$growth) and log10(lig3$growth)
t = 5.3694, df = 137.17, p-value = 3.286e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4036127 0.8742015
sample estimates:
mean of x mean of y 
 7.383824  6.744917 

t.test(log10(hum4$growth),log10(lig4$growth))

    Welch Two Sample t-test

data:  log10(hum4$growth) and log10(lig4$growth)
t = 4.6217, df = 122.92, p-value = 9.467e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.3000609 0.7496361
sample estimates:
mean of x mean of y 
 7.525765  7.000917 

t.test(log10(hum5$growth),log10(lig5$growth))

    Welch Two Sample t-test

data:  log10(hum5$growth) and log10(lig5$growth)
t = 4.9328, df = 103.98, p-value = 3.089e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.3982521 0.9337217
sample estimates:
mean of x mean of y 
 7.691692  7.025706 

hum1 = out2[out2$fire_day ==1 & out2$eco1==6,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==6,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==6,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==6,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==6,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==6,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==6,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==6,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==6,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==6,]
mean(hum1$growth)

[1] 16605316

mean(hum2$growth)

[1] 30366275

mean(hum3$growth)

[1] 44666667

mean(hum4$growth)

[1] 63374279

mean(hum5$growth)

[1] 85629090

mean(lig1$growth, na.rm=T)

[1] 2747161

mean(lig2$growth, na.rm=T)

[1] 8686182

mean(lig3$growth, na.rm=T)

[1] 15363041

mean(lig4$growth, na.rm=T)

[1] 21997987

mean(lig5$growth, na.rm=T)

[1] 26941885

t.test(log10(hum1$growth),log10(lig1$growth))

    Welch Two Sample t-test

data:  log10(hum1$growth) and log10(lig1$growth)
t = 3.1605, df = 69.967, p-value = 0.002328
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2244363 0.9921939
sample estimates:
mean of x mean of y 
 6.183231  5.574916 

t.test(log10(hum2$growth),log10(lig2$growth))

    Welch Two Sample t-test

data:  log10(hum2$growth) and log10(lig2$growth)
t = 4.8606, df = 97.58, p-value = 4.473e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4247076 1.0108416
sample estimates:
mean of x mean of y 
 7.000225  6.282451 

t.test(log10(hum3$growth),log10(lig3$growth))

    Welch Two Sample t-test

data:  log10(hum3$growth) and log10(lig3$growth)
t = 3.4215, df = 83.534, p-value = 0.0009662
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2080219 0.7855227
sample estimates:
mean of x mean of y 
 7.170545  6.673773 

t.test(log10(hum4$growth),log10(lig4$growth))

    Welch Two Sample t-test

data:  log10(hum4$growth) and log10(lig4$growth)
t = 2.9997, df = 74.706, p-value = 0.00367
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1331595 0.6597929
sample estimates:
mean of x mean of y 
 7.343272  6.946796 

t.test(log10(hum5$growth),log10(lig5$growth))

    Welch Two Sample t-test

data:  log10(hum5$growth) and log10(lig5$growth)
t = 3.1455, df = 56.335, p-value = 0.002647
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1801401 0.8117505
sample estimates:
mean of x mean of y 
 7.477756  6.981811 

hum1 = out2[out2$fire_day ==1 & out2$eco1==11,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==11,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==11,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==11,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==11,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==11,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==11,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==11,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==11,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==11,]
mean(hum1$growth)

[1] 21604386

mean(hum2$growth)

[1] 43011405

mean(hum3$growth)

[1] 71836206

mean(hum4$growth)

[1] 110978456

mean(hum5$growth)

[1] 156565228

mean(lig1$growth, na.rm=T)

[1] 12710105

mean(lig2$growth, na.rm=T)

[1] 24454241

mean(lig3$growth, na.rm=T)

[1] 39894008

mean(lig4$growth, na.rm=T)

[1] NaN

mean(lig5$growth, na.rm=T)

[1] NaN

#t.test(log10(hum1$growth),log10(lig1$growth))
#t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))

Error in t.test.default(log10(hum3$growth), log10(lig3$growth)) : 
  not enough 'y' observations

out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

mean(out1$ros_km,na.rm=T)

[1] 0.8405408

mean(out2$ros_km,na.rm=T)

[1] 1.82709

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

out1 = subset(res,cause == 1 & eco1 == 11)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11)

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ autumn northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 9))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 9))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ autumn mediterean california
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 9))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 9))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ summer northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ summer mediterean california
out1 = subset(res,cause == 1 & eco1 == 11  & (month > 5 & month<10))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11  & (month > 5 & month<10))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

mean(out1$ros_km,na.rm=T)

[1] 0.8405408

mean(out2$ros_km,na.rm=T)

[1] 1.82709

#0,0.25,0.5,1,2,3,5,7,10,20,30
hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure2_freq_v1.tif", width = 6, height = 5, units = 'in', res = 300)

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                2 

############ autumn northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 9))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 9))

hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

fr = barplot(dens, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="% fire days",ylim=c(0,80),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/S4_freq_autumn_north_v1.tif", width = 6, height = 5, units = 'in', res = 300)

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,80),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                2 

############ autumn mediterean california
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 9))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 9))

hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

fr = barplot(dens, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="% fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/S4_freq_autumn_med_v1.tif", width = 6, height = 5, units = 'in', res = 300)

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,25),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                2 

############ summer northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))

hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

fr = barplot(dens, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="% fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/S4_freq_summer_north_v1.tif", width = 6, height = 5, units = 'in', res = 300)

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                2 

############ summer mediterean california
out1 = subset(res,cause == 1 & eco1 == 11  & (month > 5 & month<10))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11  & (month > 5 & month<10))

hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

fr = barplot(dens, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="% fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/S4_freq_summer_med_v1.tif", width = 6, height = 5, units = 'in', res = 300)

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                2 

20 fastest fires

res1 = res[res$ros_km > 10 & !is.na(res$ros_km),]
length(res1$ros_km)
res1

are ROS the same for light & human under the same conditions

out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

plot(out1$vpd,log(out1$ros_km))
points(out2$vpd,log(out2$ros_km),col="red")

summary(lm(out1$vpd~log(out1$ros_km+1),na.omit=T))
summary(lm(out2$vpd~log(out2$ros_km+1)))

analysis of the first day

load data

daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/all_ignitions_V3.txt",header=T)

res=as.data.frame(daily_res)

res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$total_area =as.numeric(as.character(res$total_area))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))

res$biomass =as.numeric(as.character(res$biomass))

res = res[-1,]
res$human[res$cause ==1] =1
res$human[res$cause !=1 & res$cause !=14] =0

analysis

out1 = res[res$cause !=1 & res$cause != 14,] 
out2 = res[res$cause ==1,] 

out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1)   #1=lightning; 14=unknown; 7=arson

out1 = subset(res,eco1 == 11& res$cause !=1 & res$cause != 14)
out2 = subset(res,eco1 == 11 & res$cause ==1)

out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14 & res$mont > 5 & res$mont < 10 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1& res$mont > 5 & res$mont < 10)   #1=lightning; 14=unknown; 7=arson

t.test(out1$bi,out2$bi)
t.test(out1$erc,out2$erc)
t.test(out1$etr,out2$etr)
t.test(out1$fm100,out2$fm100)
t.test(out1$fm1000,out2$fm1000)
t.test(out1$pet,out2$pet)
t.test(out1$pr,out2$pr)
t.test(out1$rmax,out2$rmax)
t.test(out1$rmin,out2$rmin)
t.test(out1$th,out2$th)
t.test(out1$tmmn,out2$tmmn)
t.test(out1$tmmx,out2$tmmx)
t.test(out1$vpd,out2$vpd)
t.test(out1$vs,out2$vs)
t.test(out1$ws,out2$ws)
t.test(out1$biomass,out2$biomass)
t.test(out1$mean_land,out2$mean_land)
t.test(log10(out1$total_area),log10(out2$total_area))

ta = table(res$human,res$mont)
ps = barplot(ta, beside=TRUE, ylab="number of fires",xpd=T,xlab= "month", xaxt='n',ylim=c(0,300), axis.lty=1,cex.lab = 1.4,cex.axis = 1.2 )
axis(1,at=c(2,5,8,11,14,17,20,23,26,29,32,35), labels =c(1:12),xlim=c(0,36),xpd=F ,cex.lab = 1.4,cex.axis = 1.3)
legend("topright",c("human","lightning"),fill = c("grey", "black"), bty="n",cex=1.4)
text(1,285,"a)",cex=1.8)
out1 = subset(res,eco1 == 6 |eco1 == 7)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,eco1 == 11)

ta = table(out1$human,out1$mont)
ps = barplot(ta, beside=TRUE, ylab="number of fires",xpd=T, xaxt='n',xlab= "month", ylim=c(0,200), axis.lty=1,cex.lab = 1.4,cex.axis = 1.2 )
axis(1,at=c(2,5,8,11,14,17,20,23,26,29,32,35), labels =c(1:12),xlim=c(0,36),xpd=F,cex.lab = 1.4,cex.axis = 1.3 )
legend("topright",c("human","lightning"),fill = c("grey", "black"), bty="n",cex=1.4)
text(1,190,"b)",cex=1.8)

ta = table(out2$human,out2$mont)
ps = barplot(ta, beside=TRUE, ylab="number of fires",xpd=T,xaxt='n',xlab= "month", ylim=c(0,200), axis.lty=1,cex.lab = 1.4,cex.axis = 1.2 )
axis(1,at=c(2,5,8,11,14,17,20,23,26,29,32,35), labels =c(1:12),xlim=c(0,36),xpd=F,cex.lab = 1.4,cex.axis = 1.3 )
legend("topright",c("human","lightning"),fill = c("grey", "black"), bty="n",cex=1.4)
text(1,190,"c)",cex=1.8)

data_s=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/daily_mean_ros_dNBR_V6.txt",row.names=NULL)
rownames(data_s) <- c()

data_s1 = as.data.frame(data_s[,2:19])


names(data_s1) = c("lon","lat","fire","nr_day","max_land","mean_land","elevation","biomass","mean_ros","ros95","mean_dnbr","dnbr95","mean_rdnbr","rdnbr95","mean_BA_red","BA_red95","cause","size")
length(data_s1$lon)

[1] 2459

data_s1$mean_ros =as.numeric(as.character(data_s1$mean_ros))
data_s1$ros95 =as.numeric(as.character(data_s1$ros95))
data_s1$mean_dnbr =as.numeric(as.character(data_s1$mean_dnbr))
data_s1$dnbr95 =as.numeric(as.character(data_s1$dnbr95))
data_s1$mean_rdnbr =as.numeric(as.character(data_s1$mean_rdnbr))
data_s1$rdnbr95 =as.numeric(as.character(data_s1$rdnbr95))
data_s1$lon =as.numeric(as.character(data_s1$lon))
data_s1$lat =as.numeric(as.character(data_s1$lat))

data_s1$mean_land[is.na(data_s1$mean_land)]=data_s1$max_land[is.na(data_s1$mean_land)] #mean landcover gives NA when only one landcover is present
data_s1$mean_land =as.numeric(as.character(data_s1$mean_land))
data_s1$mean_BA_red =as.numeric(as.character(data_s1$mean_BA_red))
data_s1$BA_red95 =as.numeric(as.character(data_s1$BA_red95))
data_s1$biomass =as.numeric(as.character(data_s1$biomass))
data_s1$elevation =as.numeric(as.character(data_s1$elevation))
data_s1$cause =as.numeric(as.character(data_s1$cause))
data_s1$size =as.numeric(as.character(data_s1$size))

data_s1 = data_s1[data_s1$size > 50000,]
data_s1=na.omit(data_s1)

shape = shapefile("/Users/stijnhantson/Documents/data/veg_california/ca_eco_l3/ca_eco_l3.shp")
pts <- SpatialPoints(data_s1[,c("lon","lat")],P4S.latlon)
shape = spTransform(shape,P4S.latlon)
eco = over(pts, shape)
data_s1$L1CODE = eco$NA_L1CODE
data_s1$L3name = eco$US_L3NAME
data_s1$L1CODE =as.numeric(as.character(data_s1$L1CODE))

#data_s1$log_ros = log10(data_s1$mean_ros)
#data_s1$log_ros95 = log10(data_s1$ros95)
#data_s1 = data_s1[data_s1$mean_ros >0,]
data_s1$human[data_s1$cause !=1 & data_s1$cause !=14 & data_s1$cause !=17]=1
data_s1$human[data_s1$cause ==1 ]=2

data_s1$ros_km = (data_s1$ros95 *24)/1000
data_test = data_s1[data_s1$max_land == 1,]
data_test1 = data_s1[data_s1$L1CODE == 6 |data_s1$L1CODE == 7 ,]
data_test2 = data_s1[data_s1$L1CODE == 11,]

data_test1 = data_s1[data_s1$human == 1 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]
data_test2 = data_s1[data_s1$human == 2 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]

data_test1 = data_s1[data_s1$human == 1 ,]
data_test2 = data_s1[data_s1$human == 2 ,]

data_test1=na.omit(data_test1)
data_test2=na.omit(data_test2)

plot(data_s1$ros_km, data_s1$mean_BA_red,log="x",xlab="Rate-of-Spread (km/day)",ylab="Tree mortality (%)",xlim=c(0.005,30),xaxt="n",cex.axis=1.4 ,cex.lab=1.4,cex=0.8, col="black")

33 x values <= 0 omitted from logarithmic plot

axis(1,at=marks,labels=marks,cex.axis=1.4 )

points(data_test2$ros_km, data_test2$mean_BA_red,cex=0.8, col="darkgrey")
points(data_test1$ros_km, data_test1$mean_BA_red,cex=0.8, col="orange")

marks=c(0.01,0.1,1,10)
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/fig_basa_ros_all_v3.tif", width = 5, height = 5, units = 'in', res = 300)
plot(data_s1$ros_km, data_s1$mean_BA_red,log="x",xlab="Rate-of-Spread (km/day)",ylab="Tree mortality (%)",xlim=c(0.005,30),xaxt="n",cex.axis=1.4 ,cex.lab=1.4,cex=0.8, col="black")

33 x values <= 0 omitted from logarithmic plot

axis(1,at=marks,labels=marks,cex.axis=1.4 )
points(data_test2$ros_km, data_test2$mean_BA_red,cex=0.8, col="darkgrey")
points(data_test1$ros_km, data_test1$mean_BA_red,cex=0.8, col="orange")
lines(lowess(data_s1$ros_km, data_s1$mean_BA_red, f=0.41),col="black", lwd=3)
lines(lowess(data_test1$ros_km, data_test1$mean_BA_red, f=0.41),col="darkgoldenrod3", lwd=3)
lines(lowess(data_test2$ros_km, data_test2$mean_BA_red, f=0.41),col="gray40", lwd=3)
legend("topleft",legend=c("all","human","lightning"),col = c("black","darkgoldenrod3", "gray40"),lty=1, bty="n",lwd = 3, cex=1)
dev.off()

quartz_off_screen 
                2 

data_test1 = data_s1[data_s1$human == 1 & data_s1$max_land <1.5,]
data_test2 = data_s1[data_s1$human == 2 & data_s1$max_land <1.5,]
data_forest = data_s1[data_s1$max_land <1.5,]
data_test1=na.omit(data_test1)
data_test2=na.omit(data_test2)

data_s2 = data_s1[data_s1$ros_km <10,]
data_forest_10 = data_s2[data_s2$max_land <1.5,]
data_test1_10 = data_s2[data_s2$human == 1 & data_s2$max_land <1.5,]
data_test2_10 = data_s2[data_s2$human == 2 & data_s2$max_land <1.5,]
data_test1_10=na.omit(data_test1_10)
data_test2_10=na.omit(data_test2_10)

marks=c(0.01,0.1,1,10)
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/fig_basa_ros_forest_v3.tif", width = 5, height = 5, units = 'in', res = 300)

plot(data_forest$ros_km, data_forest$mean_BA_red,log="x",xlab="Rate-of-Spread (km/day)",ylab="Tree mortality (%)",xlim=c(0.005,30),xaxt="n",cex.axis=1.4 ,cex.lab=1.4,cex=0.8, col="black")

30 x values <= 0 omitted from logarithmic plot

axis(1,at=marks,labels=marks,cex.axis=1.4 )
points(data_test2$ros_km, data_test2$mean_BA_red,cex=0.8, col="darkgrey")
points(data_test1$ros_km, data_test1$mean_BA_red,cex=0.8, col="orange")
lines(lowess(data_forest_10$ros_km, data_forest_10$mean_BA_red, f=0.41),col="black", lwd=3)
lines(lowess(data_test1_10$ros_km, data_test1_10$mean_BA_red, f=0.41),col="darkgoldenrod3", lwd=3)
lines(lowess(data_test2_10$ros_km, data_test2_10$mean_BA_red, f=0.41),col="gray40", lwd=3)
legend("topleft",legend=c("all","human","lightning"),col = c("black","darkgoldenrod3", "gray40"),lty=1, bty="n",lwd = 3, cex=1)
dev.off()

quartz_off_screen 
                2 

marks=c(0.01,0.1,1,10)
tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/fig_basa_ros_forest_v4.tif", width = 5, height = 5, units = 'in', res = 300)
plot(data_forest$ros_km, data_forest$mean_BA_red,log="x",xlab="Rate-of-Spread (km/day)",ylab="Tree mortality (%)",xlim=c(0.1,15),xaxt="n",cex.axis=1.4 ,cex.lab=1.4,cex=0.8, col="black")

30 x values <= 0 omitted from logarithmic plot

axis(1,at=marks,labels=marks,cex.axis=1.4 )
points(data_test2$ros_km, data_test2$mean_BA_red,cex=0.8, col="darkgrey")
points(data_test1$ros_km, data_test1$mean_BA_red,cex=0.8, col="orange")
lines(lowess(data_forest_10$ros_km, data_forest_10$mean_BA_red, f=0.41),col="black", lwd=3)
lines(lowess(data_test1_10$ros_km, data_test1_10$mean_BA_red, f=0.41),col="darkgoldenrod3", lwd=3)
lines(lowess(data_test2_10$ros_km, data_test2_10$mean_BA_red, f=0.41),col="gray40", lwd=3)
legend("topleft",legend=c("all","human","lightning"),col = c("black","darkgoldenrod3", "gray40"),lty=1, bty="n",lwd = 3, cex=1)
dev.off()

quartz_off_screen 
                2 

plot(data_forest$ros_km, data_forest$mean_BA_red,log="x",xlab="Rate-of-Spread (km/day)",ylab="Tree mortality (%)",xlim=c(0.1,15),xaxt="n",cex.axis=1.4 ,cex.lab=1.4,cex=0.8, col="black")

30 x values <= 0 omitted from logarithmic plot

axis(1,at=marks,labels=marks,cex.axis=1.4 )

points(data_test2$ros_km, data_test2$mean_BA_red,cex=0.8, col="darkgrey")
points(data_test1$ros_km, data_test1$mean_BA_red,cex=0.8, col="orange")

lines(lowess(data_forest_10$ros_km, data_forest_10$mean_BA_red, f=0.41),col="black", lwd=3)
lines(lowess(data_test1_10$ros_km, data_test1_10$mean_BA_red, f=0.41),col="darkgoldenrod3", lwd=3)

lines(lowess(data_test2_10$ros_km, data_test2_10$mean_BA_red, f=0.41),col="gray40", lwd=3)
legend("topleft",legend=c("all","human","lightning"),col = c("black","darkgoldenrod3", "gray40"),lty=1, bty="n",lwd = 3, cex=1)

NA
NA

data_test = data_s1[data_s1$max_land == 1,]
data_test1 = data_s1[data_s1$L1CODE == 6 |data_s1$L1CODE == 7 ,]
data_test2 = data_s1[data_s1$L1CODE == 11,]

data_test1 = data_s1[data_s1$human == 1 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]
data_test2 = data_s1[data_s1$human == 2 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]

data_test1 = data_s1[data_s1$human == 1 ,]
data_test2 = data_s1[data_s1$human == 2 ,]


fast = data_s1[data_s1$ros_km > 1,]
fast_hum = fast[fast$human == 1,]

sum(fast$size)/sum(data_s1$size)
length(fast$size)/length(data_s1$size)

sum(fast_hum$size, na.rm=T)/sum(fast$size)
length(fast_hum$size)/length(fast$size)


quan = quantile(data_s1$ros_km,0.5)
fast = data_s1[data_s1$ros_km > quan,]
slow = data_s1[data_s1$ros_km < quan,]
fast_hum = fast[fast$human == 1,]

sum(fast$size)/sum(data_s1$size)
length(fast$size)/length(data_s1$size)
sum((data_s1$mean_BA_red*data_s1$size))/(sum(data_s1$size))
mean(data_s1$mean_BA_red)

sum(fast_hum$size, na.rm=T)/sum(fast$size)
length(fast_hum$size)/length(fast$size)
sum((fast$mean_BA_red*fast$size))/(sum(fast$size))
mean(fast$mean_BA_red)
sum((slow$mean_BA_red*slow$size))/(sum(slow$size))
mean(slow$mean_BA_red)

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

```